using System;
using System.Collections.Generic;
using System.Linq;

using Atomic.Libraries;
using Atomic.Libraries.Mathematics;
using Atomic.Structures.Symmetries;

namespace Atomic.Structures
{
	public static class StructureEquality
	{
		/// <summary>
		/// Compares two structures allowing small deviation in the atomic positions.
		/// This does NOT define an equivalence relation as it lacks the transitive property.
		/// The tolerance is assumed to be in real units (i.e. in Å).
		/// </summary>
		public static bool Equals(Structure structure1, Structure structure2, double tolerance)
		{
			if (tolerance < 0.0)
			{
				throw new ArgumentException();
			}

			if (structure1 == structure2)
			{
				// Avoid expensive check if they're truely equal.
				return true;
			}

			if (structure1 == null || structure2 == null)
			{
				// Case where both are null is already handled above.
				return false;
			}

			if (SpaceMatrix.Norm(structure1.Lattice.PhysicalMap - structure2.Lattice.PhysicalMap) > tolerance)
			{
				return false;
			}

			return Distance(structure1, structure2) <= tolerance;
		}

		public static bool RotationEquals(Structure structure1, Structure structure2, double tolerance)
		{
			if (Equals(structure1, structure2, tolerance))
			{
				return true;
			}

			try
			{
				SpaceMatrix rotation = structure1.Lattice.PhysicalMap * structure2.Lattice.InversePhysicalMap;
				return Equals(structure1, structure2.Rotate(rotation), tolerance);
			}
			catch (ArithmeticException)
			{
				// Not a proper rotation matrix.
				return false;
			}
		}

		/// <summary>
		/// Compares two structures allowing small deviation in the atomic positions and acting with symmetry operations.
		/// This does NOT define an equivalence relation as it lacks the transitive property.
		/// The tolerance is assumed to be in real units (i.e. in Å).
		/// </summary>
		public static bool SymmetryEquals(Structure structure1, Structure structure2, IEnumerable<SpaceGroupSymmetryOperation> symmetries, double tolerance)
		{
			if (symmetries == null)
			{
				// Only this is tested for null. The structures are allowed to be null.
				throw new ArgumentNullException();
			}

			if (tolerance < 0.0)
			{
				throw new ArgumentException();
			}

			if (structure1 == structure2)
			{
				// Avoid expensive check if they're truely equal.
				return true;
			}

			if (structure1 == null || structure2 == null)
			{
				// Case where both are null is already handled above.
				return false;
			}

			if (structure1.Sites.Count != structure2.Sites.Count || new SymmetricEquatableEnumerable<Atom>(structure1.Sites.Select(s => s.Atom)) != new SymmetricEquatableEnumerable<Atom>(structure2.Sites.Select(s => s.Atom)))
			{
				// Multiplicities of the atoms are different.
				return false;
			}

			foreach (SpaceGroupSymmetryOperation symmetry in symmetries)
			{
				if (Equals(structure1, Transform(symmetry, structure2), tolerance))
				{
					return true;
				}
			}

			return false;
		}

		/// <summary>
		/// Compares two structures allowing small deviation in the atomic positions and acting with symmetry operations.
		/// This does NOT define an equivalence relation as it lacks the transitive property.
		/// The tolerance is assumed to be in real units (i.e. in Å).
		/// </summary>
		public static bool SymmetryEquals(Structure structure1, Structure structure2, Structure lattice, double tolerance)
		{
			if (lattice == null)
			{
				throw new ArgumentNullException();
			}

			return SymmetryEquals(structure1, structure2, lattice.SpaceGroup.SymmetryOperations, tolerance);
		}

		/// <summary>
		/// Computes the maximum distance between sites in the two structures. Lattice vectors are allowed to be different.
		/// </summary>
		public static double Distance(Structure structure1, Structure structure2)
		{
			// Do mapping of sites from both structures to make sure that no site is left out from the distance.
			return 0.5 * (AsymmetricDistance(structure1, structure2) + AsymmetricDistance(structure2, structure1));
		}

		private static double AsymmetricDistance(Structure structure, Structure otherStructure)
		{
			double d0 = 0.0;

			foreach (Site otherSite in otherStructure.Sites)
			{
				Atom atom = otherSite.Atom;
				SpaceVector position = otherSite.Position.PhysicalCoordinates;

				// Find the best fitting lattice point in structure1 with the right atom.
				double d = double.PositiveInfinity;
				foreach (Site site in structure.Sites)
				{
					if (site.Atom != atom)
					{
						// Wrong atom. Just skip ahead to next site.
						continue;
					}

					// Find best possible lattice point.
					SpaceVector x = structure.Lattice.Decompose(new Position(position - site.Position.PhysicalCoordinates)).Coordinates;

					// Deviation in physical coordinates.
					d = Math.Min(d, SpaceVector.Norm(structure.Lattice.PhysicalMap * new SpaceVector(x[0] - Math.Round(x[0]), x[1] - Math.Round(x[1]), x[2] - Math.Round(x[2]))));
				}

				// Take worst fit.
				d0 = Math.Max(d0, d);
			}

			return d0;
		}

		private static Structure Transform(SpaceGroupSymmetryOperation symmetry, Structure structure)
		{
			return new Structure(structure.Basis, structure.Lattice, Transform(symmetry, structure.Sites));
			//return new Structure(new Lattice(operation.Transform(structure.Lattice.A), operation.Transform(structure.Lattice.B), operation.Transform(structure.Lattice.C)), Transform(operation, structure.Sites));
		}

		private static IEnumerable<Site> Transform(SpaceGroupSymmetryOperation symmetry, IEnumerable<Site> sites)
		{
			foreach (Site site in sites)
			{
				yield return new Site(site.Atom, symmetry.Transform(site.Position));
			}
		}
	}
}
